# Forelesning 18, fredag 10.03.2017

#1# MLE - bruk av optim             

tsim <- rweibull(50,shape=2,scale=10)
hist(tsim,breaks=10,freq=FALSE)
curve(dweibull(x,shape=2,scale=10),add=T)

# the negative log likelihood for the weibull model
#
# arguments:
#   p : a vector containing a og b
#   t : a vector containing the observed data
#
# value
#   the negative log likelihood
lnL <- function(p,t) {
  a <- p[1]
  b <- p[2]
  n <- length(t)
  -(n*log(a) - n*a*log(b) + (a-1)*sum(log(t)) - sum((t/b)^a))
}
# alternativ fuction using built in density function
lnL <- function(p, t) {
  - sum(dweibull(t,shape=p[1],scale=p[2],log=T))
}

install.packages("emdbook")
library(emdbook)
curve3d(exp(-lnL(c(x,y),t=tsim)),xlim=c(1.5,3),ylim=c(5,15),sys="persp",theta=30,phi=30,
        tick="detailed",xlab="shape parameter a",ylab="scale parameter b",zlab="likelihood")

#Hva ser det ut som estimatene skal være? 

?optim
fit <- optim(c(.1,1),lnL,t=tsim,hessian=T)
fit
fit$par


# approximate standard errors based on asymptotic theory
sqrt(diag(solve(fit$hessian)))


##Likelihood ratio tester
# Testing H0:a=1 vs H1:a<>1 based on asymptotic or approximate chisquare
## distribution of the log likelihood

## maximum log likelihood under H0 : a=1
n <- length(tsim)
-(-n*log(mean(tsim))-sum(t/mean(t)))
lnL(c(1,mean(tsim)),t=tsim) # alternativ hyp

## 2 twice the difference in the log likelihood, 2(logL1-logL0)
#NB: Disse tallene vil endre seg fra en simulering til en annen 
2*(-149-(-112))
qchisq(.95,df=1)
pchisq(2*(-140-(-104)),df=1) #Liten? Ser fra 3d plottet at a=1 gir en liten logL